The life-cycle of drift-wave turbulence driven by small scale instability 
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We demonstrate theoretically and numerically the zonal-flow/drift-wave feedback mechanism for 
the LH transition in an idealised model of plasma turbulence driven by a small scale instability. Zonal 
flows are generated by a secondary modulational instability of the modes which are directly driven 
by the primary instability. The zonal flows then suppress the small scales thereby arresting the 
energy injection into the system, a process which can be described using nonlocal wave turbulence 
theory. Finally, the arrest of the energy input results in saturation of the zonal flows at a level which 
can be estimated from the theory and the system reaches stationarity without damping of the large 
scales. 
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Designing a device capable of controlled, self-sustained 
nuclear fusion is one of the principal technological chal- 
lenges of the 21st century. Attention is now heavily fo- 
cused on the construction of the International Thermonu- 
clear Test Reactor (ITER) which aims to demonstrate the 
feasibility of controlled fusion for the first time. Of the 
myriad of engineering challenges posed by ITER, one of 
the most critical, the problem of magnetic plasma con- 
finement, relates intimately to foundational issues in tur- 
bulence theory, the topic of this Letter. One difficulty 
of magnetic plasma confinement is due to instabilities 
(such as the ion temperature gradient, drift-dissipative 
and electron temperature gradient instabilities). These 
generate waves which interact to produce wave turbu- 
lence thereby greatly enhancing transport. Heat and 
plasma then diffuse more quickly from the hot, dense 
core to the periphery resulting in a loss of confinement. 
Common features of these instabilities are that they are 
small scale, typically of the order of the ion gyro-radius or 
below, and principally generate meridional waves (with 
wave crests along the minor radius of the torus). 

Significant progress in the confinement problem came 
with the experimental discovery pQ that, under the right 
conditions, the turbulent transport coefficient suddenly 
drops. This effect, called the Low-to-High (LH) confine- 
ment transition, can be exploited to design more effi- 
cient confinement. Although traces of the LH transition 
were replicated in numerical simulations [5] , including full 
gyro-kinetic calculations [3], the underlying mechanism 
was initially poorly understood since realistic numerical 
simulations of tokamaks are universally incapable of re- 
solving the full range of scales involved. It is now gen- 
erally accepted [I] that the LH transition is triggered by 
the self-organisation of coherent zonal flows (ZF's) in the 
poloidal direction which act as transport barriers since 
they act to shear away turbulent fluctuations which try 
to cross them. Current research is heavily focused on 
understanding and parameterising these ZF's and their 
interaction with the drift wave turbulence at the experi- 



mental, numerical and theoretical levels. 

In this Letter we present a simple theoretical model of 
the confinement mechanism described above. Our model 
uses the Charney-Hasegawa-Mima (CHM) equation [5], 
the simplest single-fluid description of a strongly mag- 
netised plasma. This equation is rarely used in mod- 
ern plasma research because drastic simplifications are 
required to obtain it from the true dynamical equations. 
On the other hand, simple models often repay in terms of 
analytical tractability and conceptual clarity what they 
lose in terms of direct applicability. This is certainly the 
case here. Our model gives a semi-analytical description 
of the full life-cycle of turbulence generated by a small- 
scale primary instability. Before presenting details, we 
outline the main results. Meridional waves - drift waves 
in the case of the CHM model - are generated by a pri- 
mary instability at the scale of the ion gyro-radius. The 
subsequent evolution has three main stages. In stage 1, 
waves grow until they generate seeds of a ZF via a sec- 
ondary modulational instability 6-8 . The ZF is large- 
scale compared to the waves. In stage 2, the ZF's grows 
by direct interaction with the small scales rather than by 
a Kraichnan-style inverse cascade which would involve lo- 
cal transfer of energy between scales. This is because the 
nonlocal interaction in drift wave turbulence is intrinsic, 
i.e. it typically develops in time even for broadband ini- 
tial spectra. However, the fact that the primary instabil- 
ity produces a relatively narrow initial drift-wave spec- 
trum helps nonlocality. Due to nonlocality, the kinetic 
energy of the small scales is not conserved as it would 
be in a local cascade. In fact, the growing ZF removes 
energy from the turbulence since the combined energy 
of large and small scales is conserved. The wave turbu- 
lence which gave birth to the ZF in stage 1 is thus sup- 
pressed by its growth. In stage 3, this negative feedback 
reaches a steady state in which suppression of turbulence 
by the large-scale ZF exactly balances its generation by 
the small scale instability. We note that Eq. ([!]) also has 
a geophysical interpretation which we do not discuss (see 
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references in [5]). 

The CHM model does not contain small scale instabil- 
ities of the kind mentioned above. In what follows, we 
model these instabilities with a multiplicative linear forc- 
ing term applied to selected wavenumbers. We solve the 
following equation on a plane with the x-axis along the 
plasma density gradient Vno and the y-axis transverse 
to this: 

(d t - C) (p 2 A^ - V>) + d y i> - p 2 J[^, Atp] = 0. (1) 

Here ip(x,y,t) is the streamfunction, A is the 2D Lapla- 
cian, J[ip, Atp] = d x ipd y Aip — d y ipd x Aip and £ is a lin- 
ear operator having Fourier-space representation £k = 
7k — v m k 2m with 7k > modeling the instability and 
v m a hyperviscosity coefficient. Here, p is the ion Lar- 
mor radius at the electron temperature, = ujbP 2 /L u 
is the drift velocity, is the ion gyro- frequency and 
L n = no/ 1 Vno | is the density profile length-scale, p and 
v* are both taken to be constant. For the instability forc- 
ing, 7k, we seek the simplest form which retains the key 
features of the relevant plasma instabilities mentioned 
above, namely that they act at small scales and primar- 
ily generate meridional waves. For these purposes it suf- 
fices to take 7k = 7<5k,k;, with k, = (0,p _1 ) and 7 > 0, 
which forces a single meriodinal mode at the scale of 
the gyro-radius. In principle, a more realistic form for 
7k can be obtained by considering the linear dynamics 
of a higher level model which does contain an intrinsic 
instability such as the Hasegawa-Wakatani (HW) model 
in the case of the drift-dissipative instability 2 . The 
dimensionless parameter, x — ^ri is a measure of nonlin- 
earity of the system assuming that the linear instability 
saturates at a level determined by the amplitude of the 
nonlinear term in Eq.([T]). Eq. ([!]) has exact plane wave 
solutions ip = exp(— i(u>kt — k • x) + 7k^ — v m k 2m t) with 
dispersion relation, u>k — —v*k y /(l + p 2 k 2 ). 

Without instability or dissipation (C = 0), Eq. ([!]) 
has two invariants: E = \ J[p~ 2 ip 2 + (V-0) 2 ] dx and 
Q = !J>~ 2 (Vt/0 2 + (V 2 ?A) 2 ]dx. These are closly re- 
lated to the energy and enstrophy in the 2D Euler equa- 
tions. 2D Euler turbulence exhibits an inverse cascade 
of energy from small scales to large. Likewise, CHM tur- 
bulence transfers energy to large scales. In the weakly 
nonlinear limit, however, CHM turbulence differs from 
2D Euler turbulence in some important respects. Firstly, 
in this limit, the energy transfer turns out to be nonlocal 
in scale [9] meaning that it is dominated by interactions 
between the smallest and largest scales. Secondly, there 
is an additional statistical invariant, the zonostrophy |10j . 
which renders the transfer anisotropic [101 E] ■ The the- 
ory of nonlocal drift wave turbulence was first worked out 
in [9l [12] and a detailed summary can be found in [13] . 
In this Letter, we will present a self-consistent but qual- 
itative physical argument to arrive at the same results. 

Consider a drift wavepacket with wavenumber k mov- 
ing on a background consisting of a large-scale ZF with 
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FIG. 1: (color online) Energy contained in the forced mode 
(crosses), the zonal sector (circles) and meridional sector (tri- 
angles) as functions of time for x = 3-75 x 10 -4 . The inset 
shows the saturation level of the ZF as a function of /3. 
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FIG. 2: (color online) Snapshots of the 2D energy spectrum 
at times 11.0, 11.5 and 20 for the case x = 3.74 x 10" 4 



typical wavenumber p. The ZF need not be perfect and 
can have finite p x , but we insist that it is large scale: 
p -C k and p -C 1/p so that a the WKB approxima- 
tion is appropriate [2] . The ZF is stationary in a frame 
moving in the y-direction with the drift velocity |13j . 
As a consequence of the WKB approximation, the wave 
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packet frequency must be conserved in this frame: 



ky p k 
1 + p 2 k 2 



= k 2 Wk = v = const. (2) 



Hence, as k changes due to distortion of the wavepacket 
by the ZF, it is restricted to remain on the ID curve in 
the 2D k-space given by hereafter called a u-curve. 
A typical i>-curve is indicated by the green (dashed) line 
in Fig. [2] According to WKB theory, it is waveaction of 
the small scales, nk = ek/wk, which is conserved rather 
than the energy, e^. On each w-curve, ~ fc 2 ek from 
Eq. ([2]). Since k 2 increases as a wavepacket moves along 
a v-cuvve from the source (at the fcj,-axis) to the sink (at 
high fc's), ek must decrease to keep constant. The 
energy thus lost is transferred to the ZF. 

There are two qualitatively different regimes corre- 
sponding to weak and strong ZFs. The weak regime oc- 
curs when drift wave packets live long enough to cross 
many ZF oscillations, experiencing multiple random dis- 
tortions and executing a random walk in k-space along 
the u-curves. It can be described by a diffusion equation: 
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(3) 



where d v is the derivative along the u-curve. A detailed 
quantitative derivation of Eq. ([3| is presented in |13j . 
The diffusion coefficient, D^, is proportional to the in- 
tensity of the large-scale ZF, and can be estimated as 
-Dk ~ (k x ) 2 T, where k x ~ —k y S is the rate of change of 
k x due to the ZF shear S — d x U ~ pU, and r is the 
correlation time r = l/(p|i|) = 1 / (p\dk x 0Jk\) ~ l/(pv*)- 
In these estimates we took into account that the essential 
dynamics happens at the scales k x ~ k y ~ l/p. Thus 
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p 2 v*p p 2 v t 



(4) 



For fixed D^, Eq. ([3| indicates that rik evolves indepen- 
dently on each v-curve and is controlled by the maximum 
value of Ak = 7k — v m k 2m on the curve. If Ak is positive, 
rik will initially grow exponentially. Correspondingly, if 
Ak is negative, «k will be exponentially damped. Large 
positive Ak are obtained on curves passing close to the 
maximum of 7^. In our case this is the forcing mode, 
k = (0, ki). Negative values are found away from this 
maximum, where dissipation dominates. The growth of 
rik on a i>-curve with positive Ak cannot, however, con- 
tinue indefinitely. Diffusion along a curve transfers en- 
ergy to large scales. This in turn increases Z?k so that 
the growth due to positive Ak can begin to be tempered 
by increasingly rapid diffusion of wavepackets from the 
instability region to the dissipation region. The region of 
growing shrinks until a steady state is reached (stage 
3). Those u-curves which originally had positive growth 
reduce their growth rates to zero and the spectrum on 
them dies out. The large-scale ZF in turn saturates at a 



finite intensity determined by balancing the forcing and 
the diffusion terms on the right hand side of Eq. ^ , giv- 
ing -Dk/0 2 ~ 7max- Eq. Q gives the saturation criterion: 

^rms^rms ^ ^*7max- (^) 

where rms denotes root-mean square. Counterintuitively, 
the saturation velocity increases with v* as the flow be- 
comes more linear. Eq. ^ fails for large x when tur- 
bulence is strong and the t^-term in Eq. ([I]) can be ne- 
glected. In the strong turbulence regime, a wavepacket 
cannot cross even a single ZF oscillation since rapid dis- 
tortion by the zonal shear brings its wavenumber to the 
dissipation scale. This gives the saturation criterion 
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Comparing Eqs. ^ and the crossover between the 
two regimes occurs when 7 max ~ u*p or when U ims ~ . 
The estimate ([6]), previously put forward in [TS], is valid 
only for strong ZF's with U Tms > «*. 

We performed numerical simulations using a standard 
pseudo-spectral code at a resolution of 256 x 256 with 
m = 8, u m = 1.0 x 10" 29 , p = 0.02 and 7 = 1.2. 
The initial condition was spatial white noise of very 
low amplitude. Fig. [l] summarizes results obtained with 
P = 1.6 x 10 5 giving \ = 3.75 x 10" 4 . We divided k- 
space into a zonal sector, \k x \ > |fcj,|, and a meridional 
sector, \k x \ < \k y \ (excluding the forced mode). Fig. l] 
shows the energy of each sector and of the forced mode 
as functions of time. Initially, the forced mode grows 
exponentially, and the zonal and meridional energies re- 
main at a level set by the initial noise. This is followed 
by a sharp growth of the the zonal and the off-zonal en- 
ergies once the amplitude of the forced mode gets high 
enough to facilitate nonlinear mode coupling (stage 1). 
Next, the meridional modes, including the forced mode 
itself, are suppressed while the zonal modes retain their 
energy (stage 2). Finally a steady state is reached (stage 
3) in which the instability forcing is switched off and the 
majority of the energy is zonal. 

Fig. [2] shows the time evolution of the 2D energy spec- 
trum in k-space normalized by the maximum. The forced 
mode initially dominates but such waves are modulation- 
ally unstable [SHE]- The first indication of the nonlinear 
evolution seen in stage 1 is therefore the growth from the 
background noise of those modes whose modulational in- 
stability growth rate is greatest. In the limit of weak non- 
linearity, the modulationally unstable modes concentrate 
on the resonant curve passing through the forcing mode, 
w(k) = cj(ki) + w(k 2 ) with k = ki + k 2 for k 3 = (ki, 0) 
[S]. The maximally unstable perturbations of a merid- 
ional wave have a significant zonal component (see Fig. 
2 of (3)- These are visible, along with the resonant curve, 
in the left panel of Fig. [2j They are the seeds of the ZF's 
responsible for the LH transition. The very rapid de- 
crease in meridional turbulent energy observed in Fig. [l] 
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at stage 2 is related with the ZF-induced diffusion along 
the w-curves with a diffusion coefficient growing as the 
ZF is becoming stronger. This diffusion process is clearly 
visible in the middle frame of Fig. [2] Deviations from the 
theoretical w-curve can be attributed to absence of a large 
scale separation between the forcing and the zonal modes 
assumed in the derivation of Eq. ^ . We remark that if 
there is a nonlocal interaction with small-scale ZF's then 
diffusion proceeds along different curves, the level sets of 
zonostrophy [121 113] . One such curve, passing through 
the forced mode is shown in Fig. [2] by the blue dashed 
line. A secondary ZF branch emerges on this curve, a 
point which we hope to explore in a future publication. 

At stage 3, the zonal flow saturation estimates Eqs. ^ 
and ^ are in reasonable agreement with the numerical 
measurements in the inset of Fig. [T] The saturation val- 
ues of U Tms S rms are compared against the estimated val- 
ues over a range of i>* values. For low there is good 
agreement with the strong turbulence estimate, Eq. ([6]). 
For large there is reasonable agreement with the weak 
turbulence estimate Eq. ^ (the discrepancy is a factor 
of 2 - 3). The departure at large v* from the prediction 
t^rms'S'rms oc is likely due to the fact that the scale 
separation becomes less pronounced for large . 

To conclude, we have analyzed the zonal-flow/drift- 
wave feedback loop mechanism for the LH transition in 
an idealised model of plasma turbulence driven by a small 
scale instability. The process has 3 stages: (i) ZF's are 
generated by a modulational instability of the small-scale 
modes which are directly excited by the instability, (ii) 
the ZF's feed back on the small-scale modes and suppress 
them via a nonlocal mechanism which closely follows the 
scenario predicted in [HHH] , (iii) the energy injection into 
the system stops and the zonal energy saturates at a level 
which is in reasonable agreement with the theoretical es- 
timates, Eqs. ([5]) and ([6]), obtained from the nonlocal 
weak and strong turbulence theories respectively. 

The ZF saturation mechanism described in this Letter 
is remarkable because it is achieved entirely dynamically 
without the need to remove the small scale driving insta- 
bility or to invoke a dissipation mechanism for the large 
scales. The fact that it works transparently in a sim- 
ple model like the CHM model is encouraging for more 
realistic applications since it suggests that the mecha- 
nism is quite robust. We conclude with some physical 
estimates aiming to establish the extent to which this 
mechanism could, in principle, be observable in real toka- 
maks. The strong turbulence scaling Eq. ([6]), appropri- 
ate for U > u*, has been widely used in tokamak the- 
ory since |15j . On the other hand, the weak turbulence 
estimate ^ is new. It should be used when J7 < u*. 
Zonal flows in the DIIID and JET tokamaks have been 



reported with U ~ 5 — 35 km/sec [THl HZ] and for 
both systems ranges from 5 km/sec in the core to 100 
km/sec in the pedestal region. Both the weak and the 
strong regimes should therefore be realisable. Some idea 
of the feasibility of the weak vs strong saturation regimes 
can be obtained for the most typical instability in the 
tokamak plasma core, the Ion- Temperature Gradient in- 
stability for which 7 max — vt / VR~Lt, where R is the 
major radius and Lj> is the temperature gradient length. 
The weak/strong crossover occurs when j max /p ~ u*, i.e. 
when pp ~ L n \J RLt where L n is the typical lenght scale 
of the density gradient. For ZFs which appear inside 
pedestal regions, L n ~ Lt and pp ~ 0.1. The crossover 
condition is thus R/L n ~ 100. In large tokamaks R > 1 
m and in the pedestal region L n > 1 cm. Thus, both 
the weak and the strong regimes of the ZF generation 
can be expected. For the edge region, one should use 
another estimate for "f ma x which is based on the drift- 
dissipative instability. In this case one can show that the 
weak regime should be expected when the electron mean 
free path exceeds R. 
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